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Abstract 

After reviewing the problematic behavior of some previously suggested finite interval 
spatial operators of the symmetric Riesz type, we create a wish list leading toward a 
new spatial operator suitable to use in the space-time fractional differential equation of 
anomalous diffusion when the transport of material is strictly restricted to a bounded do- 
main. Based on recent studies of wall effects, we introduce a new definition of the spatial 
operator and illustrate its favorable characteristics. We provide two numerical methods 
to solve the modified space-time fractional differential equation and show particular re- 
sults illustrating compliance to our established list of requirements, most important to 
the conservation principle and the second law of thermodynamics. 
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1. Introduction 

In practical sense anomalous diffusion can be detected by heavy tails of the resulting 
density distribution (at a given time) and by the departure from linearly evolving mean- 
square displacement for an initially concentrated plume [ESS]- Such a behavior is 
well documented for instance for the spreading of contaminants in heterogeneous porous 
media, where shortcut pathways may be present between two points in space (causing a 
departure from Fick's law) and/orparticles can be trapped, hindered at various locations 
(causing memory effects) [j,ijEja[iL 

With a common term, the behavior may be non- 
local both in space and time [3, [lCj, . 

Continuous time random walks (CTRW) serve as a small-scale conceptual model for 
describing anomalous diffusion. Of practical interest is the evolution of the density of a 
cloud of walkers on the macroscopic scale that is ultimately determined by the statistical 
characteristics of the jump lengths and waiting times on the microscopic level. The so 
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called Levy nights are dominated by rare but large jumps and can reproduce the power- 
law tails of the spatial distributions at a given time. Allowing rare but long waiting times 
can also lead to marked departure from the scaling law of diffusion. The flux computed 
on the scale of particle motion, depends on the parameters of the random walk and on the 
density at the considered time (in the Markovian case) or on the complete density history 
(in the case of memory effects included). The passage from microscopic to macroscopic 
scale is performed b y le tting the characteristic length and time of the particle motion 



tend to zero On unbounded domain, the resulting macroscopic behavior is 

conveniently described by the space-time fractional diffusion equation [15| . The passage 
also results in the generalization of Fick's law 17 1. 



1.1. Brief summary of known results for the unbounded domain case 

The one-dimensional space-time fractional diffusion equation is written as 

g0 Qa 

-^-ttu(x, t) — a—. — ^—u(x, t) — go < x < oo, t > (1) 
dtp d\x\ a y ' w 

where a is positive constant with dimension [L a /T°]. (The dimension of u is [l/£] 
because it is understood as one dimensional density of a countable quantity.) 

The fractional time derivative is taken in the Caputo sen se I I 81] - The notation for 
the operator g ^ a was introduced by Saichev and Zaslavsky [l9|. It is understood as 

the application of the (symmetric) Riesz derivative operator with respect to the 

space variable x. The Riesz derivative is defined through the Liouville-Weyl fractional 
derivatives: 

.[D«f + D?Lf], < a < 2, a ^ 1 



2cos(7tq/2) I 

f{x) = { -£H(f;x), a=l (2) 



d\x\ a 



where D± are called the left- and right Liouville-Weyl derivatives: 

+ ~ F(m - a) dx™ J.^ (x - £)a-m+i ' 171 — \ a \ W 

D a _ (-i) m d™ r 00 fjtw 

T(m - a) dx m J x (£_ - x) a ~ m+1 ' ^ ' 

and H denotes Hilbert transform. (See a more detailed discussion, for instance, by 
Chechkin at al. [1| ) 

For the unbounded case the following Cauchy problem can be stated: solve ([T|) for 
a given parameter set {0 < a < 2, < (3 < 1, a > 0} augmented with the initial 
condition u(x,0) = fi(x), where fi is a probability density function. The solution of the 
Cauchy problem can be obtained by the Laplace-Fourier approach, probably first used 



in this context by Montroll and Weiss [2l|. The transforms serve double purpose: they 
provide a better understanding of the operators involved and also lead to the solution for 
particular cases. In fact, the /?-order Caputo derivative (0 < /3 < 1) is the generalization 
of the first derivative via Laplace transform: 

L (^ s ) =s^ 1 [sL(f;s)-f(0)} (5) 
V ' 2 



and the a order Riesz derivative (0 < a < 2) is the generalization of the second derivative 
via Fourier transform: 

KI^H = - mqf(/;w) (e) 

The fundamental solution (spreading of an initial Dirac delta) can be obtained by the 
Laplace-Fourier method and can be given in terms of well investigated special functions. 
Though various representations are available, their equivalence has been well established 
[22I [23j | . The remarkable scaling property of the fundamental solution can be stated as: 

u(x,t) = (at)-^M[-^L;^] (7) 

where M denotes the Mainardi (or M) -function given by 

1 00 (—P\ n 

««;,)^E „ !r[ _^ (1 _ ff)1 (8) 

(By including the factor 1/2 we ensure that the integral of (0 over the x-axis is unity.) 
While computability of the M-function is far from trivial, it has been basically resolved. 
For instance, Mathematica can calculate it with any desired accuracy from its definition 
above, provided the \i parameter is passed as a rational fraction. Some known special 
cases, see e.g. [24] can be reproduced symbolically with Mathematica: 

(9) 
(10) 
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where Ai stands for the Airy function and and Ai' for the Airy-prime function. 
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1.2. Anomalous spreading on a finite interval 

For Fickian diffusion the finite domain model consistent with the first law (conserva- 
tion principle) is obtained by requiring zero flux at the two endpoints of the considered 
interval. This translates to the well-known homogenous Neumann boundary conditions 
for the traditional {(3=1, a — 2} diffusion equation. Mapping the same physical require- 
ment to mathematically treatable objects for the space-time fractional partial differential 
equation has turned out to be extremely challenging. 

Numerical experimentation followed two complementary approaches: Monte-Carlo 
(Langevin) simulation of random walk and finite difference approximation to the solution 
of a space-time fractional differential equation on bounded domain. The first is easier 
to conduct and always results in physically meaningful results (including conservation, 
if reflective walls are applied), but those results are difficult to use in a practical sense. 
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The second approach has grown mature during the years. Here we will rely on the well 
documented "matrix approach" suite by Podlubny et al. [25[ that provides a general 
framework to the numerical solution of partial fractional differential equations. 

However, the problem is deeper than purely finding an adequate numerical method. 
In the presence of wall(s) of various properties the spatial operator itself needs to be 
modified, see e.g. [H, Ell- 
in this work we are looking for the solution of the Cauchy problem for the fractional 
partial differential equation 

■^gu(x, t) = a- r-j^ufo t), < x < 1, t > (12) 

when a > and the initial condition u(x,0) = fi(x), < x < 1 is a probability density 
function. The fractional time derivative of order (0 < (3 < 1) is in the Caputo sense. 
We added a subscript "mod" to the spatial operator of order (0 < a < 2), because the 
definitions (J2HU) require the extension of the u(x, t) function defined on the interval [0, 1] 
to a left function fi(x) defined on (— oo, 1] (or, strictly speaking, at least on (-co, x],) and 
to a right function f r (x) defined on [0, oo). Notice that these auxiliary functions need 
not be probability distributions. For brevity, we will call the set of choices we make in 
creating these extensions a "prescription". Various prescriptions will give rise to various 
finite domain Riesz operators and will ultimately define the characteristics of the solution 
of (HU). 

In the following section © we illustrate the problematic behavior of some previously 
suggested prescriptions and create a "wish list" leading toward a new spatial operator. 
The subsequent section ([3]) introduces the new prescription suitable to treat anomalous 
spreading on a bounded domain and describes its main characteristics. The ^ provides 
two numerical methods to solve the modified space-time fractional differential equation 
and shows particular results illustrating the key issues. We finish the paper with summary 
and conclusions. 



2. Finite domain approaches 

The mathematically straightforward prescription to create fi{x) and f r {x) is padding 
the function with zero from both sides: 



u(x,t), 0<x<l 
0, -oo < x < 



«*> = {o, (l fi,°< £ : £1 W 

With such a prescription, the operator will yield the left finite Riemann-Liouville 
fractional derivative aD^u(x,t) and the operator D°_ will yield the right finite Riemann- 
Liouville fractional derivative x Dfu(x,t), where the integrals are defined already only 
over the appropriate part of the interval [0,1], see [28]. In most numerical calculations 
published so far such a prescription has been used. It is also the default in the matrix 
approach. Starting from an fi(x) probability distribution obeying /,(0) = fi(l) = 0, it 
seems reasonable to augment the problem with the two Dirichlet boundary conditions 
it(0, t) — u(l, t) — and solve it numerically by finite differences. 
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It is well known however, that the solution with Dirichlct conditions will not satisfy 
the first law, even in the case of j3 = 1 and a — 2. Fig 1. shows the solution of a well 
documented problem with data {(3 — 1/2, a = 3/2, a = 1, fi(x) = &x(l — x)} obtained 
via the matrix approach suite. (Notice that the actual a parameter passed on should 
be a' = y/2 instead of a = 1 because in that suite the Riesz derivative is understood 
without the factor cos(7ra/2).) In the illustration we also show the evolution of the first 
integral of the spatial density that is the fraction of substance remaining in the finite 
domain. As it is obvious, material is lost during the process. Replacing the two boundary 
conditions by " fractional derivative of order a — 1 equal to zero" condition - motivated 
by some interpretation of Fick's law - does not help either. The problematic behavior 
of prescription (fTB"]) has been repeatedly discussed, for instance, in the groundwater 
literature (2J|. 

Here we show, that no boundary condition can be found to reconcile the non-physical 
nature of prescription (|13[) . To this end we introduce a small change into the problem 
specification. Instead of requiring something at the two boundaries, we require that the 
numerical solution preserve the two important characteristics of the initial distribution: 
the first integral over [0, 1] is equal to unity and it is symmetric. To satisfy these condi- 
tions instead of the two Dirichlet boundary conditions we require J Q u(x,t)dx = 1 and 
u(0,t) = u(l,t). (The second one is obviously necessary, but it turns out to be sufficient 
as well.) The two conditions are easily passed on to any finite-difference method, in 
this case to the matrix approach suite [25[. Summarized in Fig. 2 are the results for 
fi{x) = 6x(l — x). Forcing the model to satisfy the first law, we lost compliance to 
the second law. The normalized entropy Q^ILi u i In Mi/ In —, where n is the number of 
spatial mesh points) is decreasing with time. 

2.1. Turning to Caputo's idea 

Similar experiences known for practitioners have led to various ideas. For instance, 



del-Castillo-Negrete et al. [31( suggested to use a modification of the Riesz derivative 
operator, following the recipe of Caputo being so successful for the time derivative. In 
our terminology, the prescription (explicitly given here only for 1 < a < 2) 



u(x, t) - u(0, t)-x [£u(x, t)] & x=Q , 0<x<l 
0, -oo < x < 



fr{x) = f *)-„(!, *) + (!- x) [&«(*,*)] tB=1 , 0<x<l (M) 

I 0^ 1 X < ~~^. CXD 

results in the Caputo form of the finite Riesz derivative. (Notice that the prescription 
introduces a jump discontinuity between the two functions fi and f r .) Initiating the 
spreading process from the uniform distribution will leave the initial state at rest, since 
the right-hand side of (fT2")) will be identically zero. Thus the prescription resolved a 
contradiction, but now we have to face another one: starting the process from a triangular 
distribution will also leave the system at rest. 

Constructing the flux expression (Fick's law) directly with the Caputo derivative, 
Zhang et al. [321 ] introduced another variant without explicit use of the first spatial 
derivative. In our notation, the prescription will take the form: 



M*) = 



u(x,t) - u(0,t), < x < 1 
0, -oo < x < 
5 



, , _ J u(x,t) -u(l,t), < x < 1 , . 

/rW - \ 0, 1< x < oo U&j 

The advantage of prescription (|15[> is that an initial triangular distribution will not be a 
steady-state solution of problem (|12p any more. However, it is still easy to find an initial 
condition (e.g. a box function non-zero only over a part of the [0;1] interval) that will 
also be a steady-state solution - in contrast to physical intuition. Fig 3 illustrates the 
problematic behavior of prescriptions (fl3|) . (fT4|) . and l[Po]). 



2.2. Summary of desired characteristics for a finite domain model 

From the introductory numerical experiments we glean a wish list. Sought is a macro- 
scopic model in the form of fractional partial differential equation (JT2J) with the "hard" 
requirements (1-5) and "soft" ones (6-10): 

i) If the initial state is a probability distribution (non-negative and with unit area 
under the curve), this property should be preserved for any time; 

ii) If the initial probability distribution is symmetric around x = 0.5, this property 
should be preserved for any time; 

iii) The only stable steady state should be the uniform distribution; 

iv) For any non-uniform initial distribution, the entropy should monotonically increase 
with time; 

v) For integer orders, the model should reduce to known results; 

vi) Starting from a Dirac delta distribution, the solution should follow the unbounded 
fundamental solution for short times; 

vii) We still want to preserve the deep correspondence with Caputo fractional derivative 
in time and LW fractional derivative in space (allowing, however, some liberty in the 
selection of the fi and f r functions) ; 

viii) Motivation should stem from a microscopic CTRW concept; 

ix) It is desirable to have analytic solution for special cases; 

x) It is desirable to have a numerical solution method within the general framework of 



discrete fractional calculus [25 1 . 

Notice that item i) is sometimes stated as the probability preserving property, or the 
conservation principle. In this work we take the liberty to refer to it as "first law". 
Item ii) expresses the invariance with respect to directing the coordinate axis, that is we 
consider only the symmetric Riesz derivative. Items iii-iv are obviously related to the 
second law (of thermodynamics). In some applications (see e.g. stock prices, 33]) the 
listed requirements can be relaxed, but for description of spreading of material they are 
obviously necessary. 



3. The effects of walls 

3.1. A prescription by Krepysheva et al. 

Our starting point is the work of Krepysheva et al. [3J] who visualized a "reflecting" 
wall at location x — and showed that, due to its non-local character, the kernel of the 
fractional space derivative has to be modified. The rule for the hopping particle was that 
if its jump interacts with the wall, it would continue to move in the mirror direction, 
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preserving the overall length of the "initially intended" jump. In the macroscopic limit, 
the modified Riesz kernel turned out to be markedly different from the standard one 
based on finite interval left and right Riemann-Luiville derivatives. In our terminology, 



the works Krepysheva [34 . l3q ] derived the following specific prescription: Starting from 



an u(x,t) function available over the non-negative x-axis, construct: 



(x,t), < x < 1 
(—x, t), — oo < x < 



f r (x) = u(x,t), < X < oo (16) 

Calvo at al. [36[ have developed the idea further, for a rather specific geometric situation, 
when the random walk is along the perimeter of a circle. Building on these results, van 
Milligen at al. [37| recently suggested a modification of the spatial operator for our 
problem (JT2J) that involves the Hurwitz zeta function. Another extension of the idea 



based on the so-called Kolwankar-Gangal derivative - was proposed by Neel et al. [17| . 
Recently, Zoia et al. [38| also discussed the effect of walls, although not from a first law 
point of view. 

3.2. A new prescription 

This work suggests another turn in the development for the isolating two-wall case. 
We recall that the hydrodynamic limit procedure makes use of the fact that in general, 
measurements correspond to time and length scales much larger than those of particle 
motions. In our opinion, it follows that even the "extremely long" jumps of a particle 
must be shorter than the domain size. Our main idea is therefore to limit the jump 
size to one, and hence allow zero or one particle-wall interaction, but never more than 
one. This suggests a new prescription: Starting from an u(x, t) function available on 
< x < 1, construct 

MO = { z-i <£<° 

0, -oo<£<0 
«(£,*), 

/ r (0 = { u(2-£,f), l<Z<l + x (17) 
0, l + x<^<oo 

The total support of both fi and f r is always two units long but it moves with the 
location of x. The // and f r functions combined contain every function value from the 
original u(x,t) twice (one corresponding to direct jump and the other bumped from the 
wall.) Fig. 4 shows the construction of the extensions /; and f r for two specific functions 
and a specific location x = 1/3. 

3.3. The modified Riemann-Luiville- Riesz derivative 

For comparison purposes, we can cast prescription (|17|) into a more familiar form, 
using the concept of modified left and right finite-interval Riemann-Luiville derivatives: 

n a f - 1 d? " [ X frnod(Qd£ _ 

na , (-i) m rf m r +1 hq) 

x x + l,modJ p, _ , - _ x ^ a - m +l ■ y V) 
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where 



f mod (0 = /(2-on(£-3/2) + /(£)n(£-i/2) + /(-on(£ + i/2) (20) 

with the Heaviside box function defined as II (£) = 1 for |£| < 1/2 and zero otherwise. 
These finite Riemann-Luiville derivatives are based on an interval of total length 2, 
always centered at the location x where we arc interested in the derivative. Therefore, 
for non-integer a the modified finite-interval Riemann-Luiville- Riesz derivative takes the 
form 

d a 1 

~7 j i^f( X )~~n 1 TT^U-lDxmodf + xDx+l,modf) (21) 

d mo d\x\ a 2cos(7ra/2) 

Definition (|2~i"T) and prescription (|T7|) are equivalent. 

It is illuminating to compare the regularly used spatial operator (|13[) and the modified 
one (f2"Tj) for some simple functions defined over [0; 1]. Fig 5 shows the comparison for 
f{x) — x — 0.5, f(x) — (x — 0.5) 2 , and f(x) = sin(7ra;). Somewhat disappointingly, the 
modified Riesz derivative (|21[) does not eliminate singularity at the end points of the 
interval for these functions and - by and large - behaves similarly to the commonly used 
definition (|13p . However, for one function family investigated in the next sub-section, 
the difference is dramatic. 



3.4- Eigen] 'unctions and eigenvectors 

Of particular interest is the the family of functions cos(jirx), j = 0, 1, . . .. Also shown 
in Fig. 5 the comparison for f(x) = cos(37ra). We see that prescription (fT3")l leads to 
singular behavior at the endpoints - as usual. On the other hand, (|21[) is not only non- 
singular, but it almost coincides with the Fourier-Riesz derivative of the entire cos(37rx) 
function over the interval [0,1], differing only in a constant factor 1.01328.... Using 
Mathematica one can show that 

d a 

- — cos(jTrar) = c a i3 ; cos(jTnr) j = 0,l,... (22) 

U"mod \X\ 

where c a j is constant. In other words, cos(j7rx) is an eigenfunction of the operator ((2"Tj) 
with eigenvalue c aj -. Moreover, c a ,o = 0, implying that the uniform distribution has 
a modified Riesz derivative equal to zero. This property will ensure that the uniform 
distribution will be a steady state solution of (fT2"]) . We managed to obtain closed form 
expression for specific a parameters (with repeated help from Mathematica) as follows: 

C2.j = -(j*) 2 , .7 = 1,2,... (23) 
c 3 /2,, = -2(^) 3/2 C F (727), .7 = 1,2,... (24) 
Ci/2,, = -2V7^f(V27), i = l,2,... (25) 

where Cf and Sf denote the cos and sin Fresnel integrals, respectively. For other a 
values we could not obtain an explicit expression, but could still develop a simple code in 
Mathematica that calculates the eigenvalue with any required number of digit accuracy 
for an a given as a rational fraction (see Appendix). 
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4. The solution of the space-time fractional differential equation on the in- 
terval [0,1] 

We introduce two approaches. The first one uses Laplace transform in time and 
collocation in space. The second one is based on finite differences. 

4-1. The Laplace Transform - Collocation method: an example 

We illustrate this method on a simple example: {/? = 1/2, a — 3/2, a = 1, fi(x) = 
1 — cos(27rx)}. Using 6 collocation points: x c = {0, -g, |, |, |, 1} we seek the Laplace 
transform of the solution at an arbitrary point, x. Introducing the vector notation 

v(s) = {cq(s),c 1 (s),c 2 (s),c 3 (s),c 4 (s),c 5 {s)} (26) 

g(x) = {1, cos(7ra;), cos(27ra;), cos(37ra;), cos(47rx), cos(57ra;)} (27) 
e = {0,-27r 3 / 2 C F (v / 2),-4v / 27r 3 / 2 C F (2),~6v / 37r 3 / 2 C* F (v / 6), 

-167r 3 / 2 C i r(2 V / 2), -10V57r 3 / 2 C* F (%/lO)} (28) 

we represent the Laplace space solution at x as 

U{x,s)=v{s).g{x) (29) 

Then its modified Ricsz derivative of order 3/2 takes the form 

Q3/2 



drnod\x\ 3/ '- 



-U(x,s)=v(s).[eg(x)] (30) 



(with component by component multiplication between e and g) and its Caputo deriva- 
tive of order 1/2 is written as 

L (w/2 u ( x > *)> s ) = 5-1/2 t s ( v ( s )-sW) - /<(*)] ( 31 ) 

Writing the partial differential equation 

s- 1 ' 2 [s (v(s).g(x)) - fi(x)} - v(s). [eg(x)} = (32) 

at the 6 collocation points x c , we obtain a system of linear equations in the 6 unknown 
variables, v(s). The solution turns out to be 

v(s) = {-,0,- _ ,0,0,0} (33) 
s 4v / 27r 3 / 2 C F (2)\/s + s 

and hence we obtain 

r w x 1 cos(27rx) .„,. 

U(x, s) = = — — i '- (34) 

y J s 4 v / 2W 2 C F (2)yi+s V ' 

This can be inverted on s resulting in 

u{x, t) = l + cxp (327r 3 C F (2) 2 i) erf UV2tt 3/2 C f (2)V?) - 1 cos(2tt2;) (35) 

Increasing the number of collocation points does not change the solution ([55")) , it is already 
exact. 
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4-2. The fundamental solution for (3 = 1/2, a = 3/2 

We can repeat the above procedure for fi (x) = S(x— \) = 1 — 2cos(27ra;) + 2cos(47ra;) — 
2 cos(67ra;) + . . . and obtain the fundamental solution: 



u(x, t) = 1 + 2 exp (327r 3 C F (2) 2 t) cos(2ttx) ^erf [4V2TT 3/2 C F {2)Vij - 1 
2exp (2567r 3 C f (2V2) 2 ^ cos(47rir) erf (l67r 3/2 C F (2\/2)%/t) - 1 
2exp(8647r 3 C F (2\/3) 2 ^ cos(67ra;) erf (i2\Z6tt 3/2 Cf(2\/3) Vtj 



(36) 



Looking back at our wish list, we would like to check item vi) requiring the correspondence 
of the unbounded and bounded solutions at early times. The analytic solution of the 
problem stated on unbounded domain would take the form: 

Uoo (x,t) = HM(V|i- 1/3 ;l/3) (37) 

where M(£; 1/3) is given by (flC)]) . Fig. 6 compares, at an early time, t = 1CP 3 , the new 
fundamental solution (|36[) computed with 50 terms and the well-know unbounded solution 
(f3"T|) . We see that the effect of the boundary has just started to show. (Unfortunately, the 
numerical evaluation of the fundamental solution is not trivial - even with Mathematica 
- because of the extremely large exponents and hence we could not increase the number 
of terms to get rid of the small oscillations of the curve.) 

Pursuing the "analytic" solution further has other drawbacks too. The fundamental 
solution would contain Mittag-Leffler functions (in addition to the eigenvalues, c a j and 
eigenfunctions cos(j 7r x) ) in the case of a general a and (3. Moreover, any other (non 
Dirac delta) initial condition would necessitate further numerical convolution. 

One can, however, easily construct the system of linear equations for a given set 
of a, (3, a and fi(.) at any selected value of the Laplace variable s, solve the system 
numerically by Gauss elimination and substitute v(s) into (|29[) . Therefore, we have a 
way to calculate (at any specified x) the Laplace transform of the solution numerically, 
and hence we can use a numerical inversion technique [39j |. The procedure is robust, if 
care is taken to do the Gauss elimination with multiple precision - large enough with 
respect to the number of terms in the GWR algorithm [40j. We will call the method 
LT-collocation with numerical inversion. A Mathematica realization of the algorithm is 
provided in the Appendix. 

We emphasize that the procedure does not require any boundary conditions, rather 
x = and x = 1 are included in the set of collocation points. (The physical "boundary 
conditions" are taken care of within the spatial operator itself.) 

Shown in Fig. 7 is the summary of the results of the procedure for {(3 = 1/2, a = 
3/2, a = 1, fi{x) = S(x — |)}, using 51 collocation points. Since the LT-collocation 
method with numerical inversion cannot be started from the Dirac delta "function", we 
pass on the solution of the unbounded problem ([3"T]) at a very early time, t = 0.00001 
as a new "initial" condition. Then we do the numerical inversion for t — 0.00001 where 
t is the time we are interested in. The "fraction of substance still in the domain" is not 
shown in Fig. 7, because at the " initial" state it is unity (for all practical purposes) and 
hence it remains unity during all times. 
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Next we compare traditional {/3 = 1, a = 2, a = 1} and space-time fractional 
spreading {/3 = 1/2, a = 3/2, a = 1} starting from a non-symmetric initial distribution 
fi(x) — 12x(l — a;) 2 , as illustrated by Figs 8 and 9. The number of collocation points 
is kept at 51. We find that the spreading is initially faster for the space-time fractional 
case, but at late time traditional dispersion becomes faster. There is also a remarkable 
difference in the "overshoot" of the density at x = - at least for the studied time points. 
(We note, however, that such comparisons are of limited value, because the parameter a 
is not the same in the two cases, in spite of looking the same.) 

4-3. Solution by the method of finite differences 

While the LT-collocation with numerical inversion works effectively, it is still illumi- 
nating to solve the problem within the framework of discrete fractional calculus. The 
form of the modified Riesz operator (f2"Tj) suggests that a relatively small modification to 
the established matrix approach will suffice. Indeed, one has to use prescription (fTT]) to 
pad the list of available function values and add a small correction to assure that for a 
constant function the modified Riesz derivatives yield zero. 

We introduced this modification into the matrix approach suite (see Appendix for 
some details). When working with the suite, we do not use the concept of " mathematical 
boundary conditions" at all, rather we write the equations for the endpoints x = and 
x = 1 as well. Therefore, the total number of equations remains the same as the number 
of unknowns. Shown in Fig. 10 is the summary of results for our previous example {[3 = 
1/2, a — 3/2, a — 1} (that is a' — y/2) when the initial condition is fi(x) — 12ik(1 — a;) 2 
and the finite difference step sizes are Ax = 0.02 and At = 0.01. (Notice that the 
number of unknowns in the matrix approach was 5000 and the system matrix had 2.5 
million elements, hence we were limited by computer memory.) Regarding accuracy, the 
results are still behind the ones depicted in Fig. 9, but the overall correspondence is 
remarkable. In particular, the fraction of substance still in the domain has less than 2 % 
error. Not only the modified Riesz derivative is "probability preserving" but also is our 
finite difference representation of it. 

5. Summary and conclusions 

We have introduced a new version of the finite Riesz derivative. The new spatial op- 
erator - combined with the Caputo derivative in time - results in a space-time fractional 
differential equation that is well suited to describe anomalous spreading of substance in 
a finite domain. The space-time fractional differential equation satisfies our postulated 
requirements. If the initial state is a probability distribution (non-negative and with unit 
area under the curve), this property is preserved for any later time; if the initial proba- 
bility density is symmetric around x = 0.5, this property is also preserved. We could not 
prove rigorously that the only stable steady state is the uniform distribution, but all spe- 
cific analytic formulaes and numerical examples indicated so. In all our examples, starting 
from a non-uniform initial distribution, the entropy monotonically increased with time. 
For integer orders, the model reproduces the known results of Fickian-Markovian diffu- 
sion over a finite domain. Starting from a Dirac delta distribution, the solution follows 
the unbounded fundamental solution of the space-time fractional differential equation for 
short times. We could preserve the deep correspondence with Liouville-Weyl fractional 



11 



derivative in space by creating the appropriate prescription, however we are not sure 
that this is the only (or even the best) way to do it. The new approach arose from a 
microscopic CTRW concept and we could manage to provide analytic solution for spe- 
cial cases. While our preferred solution method is the LT-collocation with numerical 
inversion, we could also extend the matrix approach and hence fit our operator into the 
mainstream framework of discrete fractional calculus. In this work we focused on "pure" 
diffusion but we do not envisage any difficulty in considering simultaneous advection or 
other external potential field. 
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6. Appendix 

Calculations were done in Mathematica. Here we show some code snippets and 
results in order to ease reproduction of our results. 

The code snippet for the LT-collocation method with numerical inversion is the following: 
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vari[s_] = Table[c[j][s], {j, 0, nm}]; 

cosi[x_] = Table[Cos[(j7r) x], {j, 0, nm}]; 

symi[x_] = Table[scos[a][j][x], {j, 0, nm}]; 

equl[s_][x_] = s A ((3 — 1) (s vari[s].cosi[x] — fi[x]) — a vari[s].symi[x]; 

xc = 7V[Range[0, 1, 1/nm], maxprec]; 

{mb[s_], As[s_]} = Map[Normal, CoefficientArrays[equl[s] [xc], vari[s]]]; 

sols[s_] := sols[s] = LinearSolve[As[s], — mb[s]]; 

Uxs[x_][s_] := sols[s].cosi[a;]; 

uxt[M_][x_][t_] := GWR[Uxs[x],t,M]; 



The user has to provide (3, a, a, the function /i[x_] in addition to the integer nm (the 
number of collocation points minus one.) The number of terms in the GWR algorithm 
M (40j and the maximum used precision maxprec (we used 200) are also required. Once 
defined, the uxt[M][x][t] function can be used to calculate the solution at a specific x and t. 
The scos[a] [j] [x] expression should evaluate to — 2j 3 ^ 2 n 3 / 2 FresnelC[\/2j]Cos[jTTx} when 
a = 3/2. In general, it will be the product of an eigenvalue cei[a][j] and the appropriate 
Cos[jTTx]. While no general formula is currently available for the eigenvalue, we can 
calculate it by the following code snippet for a given rational fraction a and positive 
integer j: 



eug = EulerGamma; 

h [j]l x -} = Cos[jirx}; 

ce\[a]\j] = iV[symRmod[/i[j], eug, a]/h[j][eug], maxprec]; 

where 

symRmod[/_,x_,a_] := - — — ]■ — — (IRmod[/, x, a) + rRmod[/, x, a}); 

2Cos[an/2\ 

In the above code the left modified Rieman-Luiville-Riesz derivative is calculated from 



IRmod[f_, t_, a_] := Module[{fext, m = Ceiling[a]}, 

fext[x_] = f[x] UnitBox[a; - 1/2] + f[2 - x] UnitBox[x - 3/2] + f[-x] UnitBox[x + 1/2]; 
lf[lntegerQ[a], D[f[t], {t, a}], 

Gammalm - a ] ^[( lntegrate [(eug^rSfT^) ' {r ' eug ~ lj eug}])/ ' eug ~* *' {< ' m}]]] 
and the right modified finite-interval Rieman-Luiville-Riesz derivative is calculated from 
rRmod[f_, t_, a_] := Module[{fext, m = Ceiling[a]}, 

fext[x_] = f[x] UnitBox[a; - 1/2] + f[2 - x] UnitBox[a; - 3/2] + f[-x] UnitBox[x + 1/2]; 
If [lntegerQ[a], D[f[t},{t,a}}, 

Gamma'tm- a] J [( lntegrate f (r-eu^ll+i-m) ■ i r > eug ' eug + !}])/- eu g - *> {*. m M 

(The extensive use of "EulerGamma" is somewhat arbitrary but proved useful in the 
context of the current version - v7.0 - of Mathematica [411].) 
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Without going into details of the derivation, here we illustrate the concept of ex- 
tending the matrix approach suite of Podlubny et al. 25]. For 1 < a < 2, the current 
symmetric Riesz function of the suit (called ransym) would hypothetically provide the 
following array SR of symmetric Riesz derivatives at points {0, lh, 2h, Zh, 4/i} where the 
function values {yo, yi,U2, U3, Da} are known: 

SR Q = 2~ 1+2a (l ~a)y + 

2- 1+2a yi 

SRi = i4- 2+Q (24+12(-l + a)a)?;o- 

4 Q CM/l + 

U~ 2 + a {24+12(-l + a)a)y 2 - 
U- 1+a (-2 + a)(-l + a)ay 3 + 
|4- 2+a (-3 + a) (-2 + a)(-l + a)ay 4 
SR 2 = U- 1+a {2-a){-l + a)ay + 
|4-i+a(6 + 3(_i + a)a)|/i- 
4 a cn/ 2 + 

H- 1+a (6 + 3(-l + a)a)y 3 + 
| 4 -i+a(2 _ a )(_l + a )ay 4 
SR 3 = U- 2+a a(-6 + lla-6a 2 + a 3 )y - 
|4- 1+Q a(2 - 3a + a 2 )yi + 
4- 1+a (2-a + a 2 )y 2 - 
4 Q ay 3 + 

A- 1+a (2-a + a 2 )y 4 
SR 4 = 2- 1+2a y 3 + 

2- 1+2q (1 - a)y 4 

A symmetric Riesz function modified according to prescription (JTTJ) would hypothet- 
ically provide the following array SRM for the same input: 
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a 2 )y 3 - 



SRM = -2 2a ay + 

2- 1+2a (2 -a + a 2 ) yi - 

±2- 1+2a a(2 - 3a + a 2 )y 2 + 
| 2 -3+2 Q( -_3 + a j a ( 2 _ 3q, + u , „ 

SRM = 4- 1+a (2-a + a 2 )y Q - 

U- 1+a a(lA-ia + a 2 ) Vl + 
U-2+ a / 24 _ i 8a + 23a 2 - 6a 3 + a 4 )y 2 + 
U- 2+a { - 48 + 92a - 58a 2 + 16a 3 - 2a 4 )j/ 3 + 
§4- 2+Q ( - 6a + 11a 2 - 6a 3 + a 4 )y 4 

SRM = -i4- 1+Q a(2-3a + a 2 )y + 

|4- 2+Q ( - 48 + 52a - 70a 2 + 20a 3 - 2a 4 )j/ 2 + 
f 4 -2+^ 24 _ lga + 23a 2 _ 6a 3 + a 4) y3+ 

|4- 2+Q ( - 8a + 12a 2 - 4a 3 ) y 4 
SRM = |4- 2+Q a( - 6 + 11a- 6a 2 +a 3 )y - 
U- 3+2a (24 - 46a + 29a 2 - 8a 3 + a 



S*i?M 



1 — \^zi — ioa -f zya — oa -p- a )yi~t 
| 4 -2+ Q ^ 24 _ 18q + 23a 2 _ g a 3 + a 4}y 2 + 

|4- 2+Q ( - 56a + 12a 2 - 4a 3 )j/ 3 + 
|4- 2+Q (24- 12a + 12a 2 )y 4 
= |2- 3+2a ( - 24 + 50a - 35a 2 + 10a 3 - a 4 )y + 
f 2~ 3+2tt a( - 6 + 11a - 6a 2 + a 3 )yi 
-i2- 1+2tt a(2-3a + a 2 )y 2 + 
2- 1+2Q (2-a + a 2 )y 3 
-2 2a ay A 

One can easily check by substitution, that when y — y 1 = y 2 = ys = y 4 , each 
derivative is zero SRM = SRMi = SRM 2 = SRM 3 = SRM 4 = 0, and when y = 
1,2/1 = 0.5, j/2 = 0,?/3 = -0.5,2/4 = -1, the sum is zero SRM + SRM X + SRM 2 + 
SRM 3 + SRM 4 = 0. 



7. Figure captions 
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Figure 1: Cauchy problem with Dirichlet boundary conditions. Results from matrix method {/3 = 
1/2, a = 3/, a = 1, a' = \fi). Upper triangles: initial condition, fi(x) = 6x(l — x), < x < 1. 
Lower triangles: solution at time t = 1. The solid lines show density distributions at intermediate times 
{0.01,0.02,0.05,0.1,0.2,0.5}. The inserts show the fraction of substance still in the domain and the 
evolution of normalized entropy (see Appendix). 




Figure 2: Cauchy problem with fixed amount of substance and symmetry condition. Results from 
Matrix method {/3 = 1/2, a = 3/2, a = 1, a' = \/2} . Upper triangles: initial condition, fi(x) = 
6x (1 — x), < x < 1. Lower triangles: solution at time t = 1. The solid lines show density distributions 
at intermediate times {0.01,0.02,0.05,0.1,0.2,0.5}. The inserts show the fraction of substance still in 
the domain and the evolution of normalized entropy. 
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Figure 3: a) Initial condition uniform distribution (solid), b) triangular distribution (dashed), c) narrow 
uniform distribution (dotted). Prescription (13) will either contradict first or second law for a). Pre- 
scription (14) will contradict second law for b) and c). Prescription (15) will contradict second law for 
c). 



f(x) = .r-0.5. Osx: 




f(x) = [,-l).v. tol 



Figure 4: Illustration to prescription (17). Construction of /; and f r (dashed) from f(x) given over 
< x < 1 (solid) for a specified x = 1/3 (dotted). Dashed line, left from x = 1/3 corresponds to 
right from x = 1/3 to f r . 
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Figure 5: Modified finite Riesz derivative - prescription (17) (solid) and standard Riesz derivative 
- prescription (13) (dashed) for various f(x), < x < 1 functions when a = 3/2. For the cosine 
function we also show the Fourier-Riesz derivative of the unbounded domain function — Equation (6) 
(dotted); it virtually coincides with the modified Riesz derivative - prescription (17), apart from a factor 
of 1.01328 . . ., see Equation (24). 
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Figure 7: Fundamental solution of the bounded domain problem {/3 = 1/2, a = 3/2, a = 1} calculated 
from the LT-collocation method with numerical inversion. Upper triangle: initial condition, taken as 
Equation (37) at a very early time t = 0.00001. Lower triangle, solution at time t = 1. The solid lines 
show density distributions at intermediate times {0.01, 0.02, 0.05, 0.1, 0.2, 0.5}. One insert shows the 3D 
surface and the other the evolution of normalized entropy. 
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Figure 8: Results from LT-collocation method with numerical inversion {/3 = 1, a = 2, a = 1} coincide 
with known results. Upper triangles: initial condition, and fi(x) = 12x(l — x) 2 , < X < 1. Lower 
triangles: solution at time t = 1. The solid lines show density distributions at intermediate times 
{0.01, 0.02, 0.05, 0.1, 0.2, 0.5}. The inserts show the 3D surface and the evolution of normalized entropy. 




Figure 9: Results from LT-collocation method with numerical inversion for {/3 = 1/2, a = 3/2, a = 1}. 
Upper triangles: initial condition fi(x) = 12x(l — x) 2 , < X < 1. Lower triangles: solution at time 
t = 1. The solid lines show density distributions at intermediate times {0.01,0.02,0.05,0.1,0.2,0.5}. 
The inserts show the 3D surface and the evolution of normalized entropy. In contrast to Fig. 8, there is 
less "overshoot" at x = and the entropy values are higher, except for t = 1. 
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Figure 10: Results from (extended) the matrix approach using the modified discrete operator (see 
Appendix) without any specific boundary conditions for {/3 = 1/2, a = 3/2, a = 1, a' = \fi}. Upper 
triangles: initial condition fi(x) = 12x(l — x) 2 , < x < 1. Lower triangles: solution at time t = 1. The 
solid lines show density distributions at intermediate times {0.01,0.02,0.05,0.1,0.2,0.5}. The inserts 
show the fraction of substance still in the domain and the evolution of normalized entropy. Notice 
the slight deviation from Fig. 9, because of the limited accuracy of the finite difference approach with 
Ax = 0.02 and At = 0.01. 
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